Excited state calculations in solids by auxiliary-field quantum Monte Carlo 
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We present an approach for ab initio many-body calculations of excited states in solids. Using 
auxiliary-field quantum Monte Carlo, we introduce an orthogonalization constraint to prevent col- 
lapse of the stochastic Slater determinants in the imaginary-time propagation. Trial wave functions 
from density-functional calculations are used for the constraints. Detailed band structures calcu- 
lated for standard semiconductors are in good agreement with GW and experimental results. For 
the challenging ZnO wurtzite structure, we obtain a fundamental band gap of 3.26(16) eV, consistent 
with experiments. 
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The accurate calculation of excited states in extended 
systems is a leading challenge in modern electronic struc- 
ture theory. Density functional theory (DFT), which usu- 
ally gives good results for ground state properties of ma- 
terials without strong electron correlation, suffers from 
the well-known band-gap problem even for simple semi- 
conductors [H-Q. The use of hybrid functionals, which 
incorporate a portion of exact exchange from Hartrce- 
Fock (HF) theory, has led to significant improvements for 
semiconductors with small to medium gaps, but not for 
large gap systems Time-dependent DFT and many- 
body perturbative approaches have shown considerable 
promise [5j. The GW method is perhaps the simplest 
and most widely used of the latter and has led to dra- 
matically improved results, largely for simple sp bonded 
materials Q. It is less successful in materials that are 
more strongly correlated. In wurtzite structure ZnO, for 
example, the accuracy of GW quasi-particle excitation 
energies is still a matter of some controversy, involving 
many factors such as the choice of pseudopotential, ap- 
proximate exchange-correlation functional, and choice of 
plasmon-pole model • A general approach which al- 
lows accurate calculations of electronic excitations across 
a wide variety of solids is thus very much in need. 
Quantum Monte Carlo (QMC) [Hi ~ 



11] 



is a non- 



(QMC) 

perturbative, many-body computational method which 
is uniquely capable of scaling up to large system sizes 
and which in principle does not involve empirical param- 
eters. Ground-state calculations with several flavors of 
QMC have seen significant recent development lfjl- 12\. 
Although some excited states calculations have also been 
performed, for example, calculations in molecules 13l.ll4| 
and diffusion Monte Carlo (DMC) calculations of excita- 
tion energies of silicon and diamond at some high sym- 
metry fc-points HI HI, the capability of QMC to treat 
excited states in general is much less developed. A major 
reason for this is the intrinsic difficulty of maintaining or- 
thogonality with the lower-lying states when the targeted 
many-body excited state is being represented stochasti- 
cally in an imaginary-time projection method. 

The auxiliary-field quantum Monte Carlo (AFQMC) 
method [H[ offers a new framework for addressing this 



difficulty and doing excited state calculations in solids. 
The random walks in AFQMC take place in a space of 
Slater determinants. The single-particle orbitals in each 
Slater determinant are expressed explicitly in terms of a 
chosen one-particle basis. Thus, in addition to orthogo- 
nalizing the single-particle orbitals with each other (as is 
done in ground-state calculations), one can orthogonalize 
them with empty (hole) orbitals to remove contamina- 
tion in the projection. For example, as an approximate 
constraint, the occupied orbitals can be orthogonalized 
with unoccupied orbitals defined by a trial wave func- 
tion. Recent developments in ground state calculations 
have demonstrated that the AFQMC framework, with 
properly formulated sign or phase constraints using trial 
wave functions, allows sign-problem- free simulations with 
high accuracy across a wide range of both strongly corre- 
lated model Hamiltonians 17| and real material systems 
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In this paper, we show how the AFQMC approach 
can be formulated to accurately calculate excited states 
in extended systems. The ground-state method is aug- 
mented by an orthogonalization constraint with virtual 
orbitals to prevent the collapse to lower-lying states in the 
imaginary-time projection. Simple trial wave functions 
from HF or DFT are used for the phaseless and orthogo- 
nalization constraints. We illustrate the method by cal- 
culating the detailed band structures in diamond struc- 
ture silicon and carbon, which are, respectively, small 
and large band gap semiconductors. We then apply the 
method to determine the bandgap in ZnO, which has 
drawn much recent attention and which presents a sig- 
nificant challenge, with high-lying strongly correlated 3d- 
bands interacting with 3s and 3p semicore states. 

In AFQMC, we target an excited state of the many- 
body Hamiltonian H with the imaginary-time projection 

1**09)) K e ~ fiH IV't)' wnere is a tria l wave func- 

tion of the excited state, and j3 = NAt with N the 
simulation time step and At the Trotter step size. The 
wave function at any imaginary-time can be thought of as 
|**08)) = Y^d, c 4> \<h) ® \4>ir)- Each Slater determinant in 
the sum is a random walker, and their distribution gives a 
statistical representation of the coefficients c^. Explicitly, 



\(pf) has the form (fa,--- ,fa,--- , fan): where each fa is 
an orbital (expressed in terms of the chosen single-partle 
basis) that evolves with ft, and m denotes the number of 
spin-f electrons. (Below we will suppress the spin index 
where possible; the spin-^ component has a similar form.) 
Thus the AFQMC process resembles the propagation of a 
population of mean-field solutions in time-dependent ex- 
ternal fields. Periodically in the propagation the orbitals 
within each random walker are orthogonalized with each 
other to ensure that the signal for fermionic antisym- 
metry is not lost in numerical noise. This step is also 
needed in ground-state calculations, and indeed even in 
mean-field calculations. The excited state energy E* can 
be calculated by the mixed estimator (l4| : 



E*(ft) 



(^\H\V*(ft)) 
<^|**G8)> 



(1) 



which converges to the exact result at ft — > oo. If \tp^} 
belongs to an irreducible representation of the symme- 
try group of H different from that of the ground state, 
Eq. (jTJ) is exact for the lowest excited state of that sym- 
metry. Otherwise, imaginary-time projection is consider- 
ably more challenging, since the propagation will tend to 
collapse to the ground state (or other lower-lying states) 

The AFQMC formalism, however, provides the abil- 
ity to prevent the collapse by imposing an additional 
orthogonality constraint naturally, using virtual or- 
bitals. For a concrete illustration, let us consider a 
targeted many-body state corresponding to the "sin- 
gle" excitation of replacing the i 
a conduction orbital labeled by 
Each random walker is still an m-electron Slater de- 
terminant. For the purpose of orthogonalization, 
however, we will regard it as the extended ordered 

list (01,- ■• , fa-i, fa, , • • • ,(j>mAm+l,--- ,fa-l,fa)- 

Any orbital denoted by '-' is a virtual orbital whose only 
role is in the orthogonalization of the m occupied orbitals. 
As an approximation, we will use the corresponding or- 



th valence orbital by 
j (thus j > m). 



bitals in |i/>t), 



i.e. 



The choice of \i/jt) is further 



discussed below. A Gram-Schmidt orthogonalization on 
this extended set ensures that the following orthogonal- 



ity conditions are obeyed: i) (fa\faj) = 0; ii) 





for k e [i + l,m]; hi) (fa\fa) = for k £ [m + l,j-i\. If 
the valence state fa is degenerate, its partners are not in- 
cluded in the constraint, which is consistent with ground 
state AFQMC for an open-shell system. Similarly, any 
degenerate partners of the conduction state fa are not 
included. After the Gram-Schmidt step, only the m oc- 
cupied orbitals are retained and included in the propaga- 
tion and measurement, until the next time for orthogo- 
nalization when the fas are re-inserted and the procedure 
repeated. The above procedure generalizes straightfor- 
wardly if the targeted excited state corresponds to "dou- 
ble" excitations or beyond. 




FIG. 1. (Color online) Illustration of AFQMC projections 
of the ground- and excited-state (i = 2 — ^ j = 5) energies 
[Eq. 0] in silicon at k = (0.3,0,0). Phaseless AFQMC re- 
sult for the ground state (red diamonds) is compared to exact 
free-projection (black circles). For the excited state, free- 
projection (blue squares) collapses. Result from the new 
method with phaseless and virtual orthogonalization con- 
straint is shown in red triangles, and compared to the GW 
excitation, indicated by violet dashed lines [HJ . 



We take simple trial wave functions directly from DFT 
calculations for the phaseless constraint ll| and for the 
orthogonalization. The starting point for constructing 
IV't(^)) & t the selected k point in the Brillouin zone (BZ) 
is the single Slater determinant DFT wave function for 
the corresponding ground state, \ipT(k))- The orbitals 
fay i are given by the eigenfunctions at k based on a well- 
converged density integrated over fc-points. For a given 
excitation of spin-er from an occupied orbital i to an un- 
occupied orbital j, we replace fa[ a by faj a - For insu- 
lators, the ground states are closed shell configurations, 
so this type of excited state Slater determinant will be 
spin-contaminated in general. To avoid this, we form a 
two-determinant singlet wave function: 



IV4) = ( a ]^ a it + ajj.^lV'T) 



(2) 



where <zj and a 



j and di are creation and destruction opera- 
tors for unoccupied (conduction) and occupied (valence) 
states, respectively, and we assume that the spatial part 
of the valence and conduction orbitals are spin indepen- 
dent in \iPt)- In degenerate cases, the trial wave func- 
tion is constructed by considering all possible promo- 
tions among these orbitals, and the coefficients of each 
determinant are set equal. All our calculations are per- 
formed in the primitive cell. In supercells, the folding 
of bands creates additional mixing of crystal momentum 
eigenstates, whose decoupling will need further study. 

Figure Q] illustrates our method in fee silicon. The 
excited state corresponds to the excitation from band 
i = 2 to the lowest lying conduction band j = 5 at 
k = (0.3, 0, 0). The trial wave function \ipr) was obtained 
from DFT with the local-density approximation (LDA), 



while for excited state \tpV) is the singlet trial wave func- 
tion of Eq. ([2]) . For comparison, AFQMC results are also 
shown from free-projection, which does not impose the 
phaseless constraint [llj and is exact for the ground state, 
but is eventually overwhelmed by the fermion sign prob- 
lem. (Large runs with ~ 2 x 10 6 walkers were used in the 
free-projection calculations.) The ground state phaseless 
AFQMC results are in excellent agreement with the exact 
calculation, indicating very small systematic error from 
the phaseless approximation, consistent with earlier re- 
sults [22j . In contrast with the ground state, exactness 
is not ensured in free projection for the excited state; a 
more severe sign problem and onset of collapse to the 
ground state is seen. The phaseless and the additional 
orthogonality constraints stabilize the excited state cal- 
culation and yield an accurate excitation energy, which 
is within ~ 0.3 eV of the GW result 21] after correction 
of finite-size effects, as we discuss next. 

Excitation energies from valence state i to conduction 
state j at any selected BZ k point are calculated as the 
difference between the AFQMC total energy and the cor- 
responding ground-state energy 
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(3) 



Because the many-body calculations are performed in 
finite-size (FS) simulation cells, the energies have FS er- 
rors which must be removed lj| 23|-26[. FS corrections 
to the thermodynamic limit of an infinite-sized supercell 
can be obtained as post-processing corrections [3, [HI 
from lower-levels of theory for both one-body (lb) and 
two-body (2b) effects. In the present excited state cal- 
culations, however, the creation of the electron-hole (eh) 
pair results in an additional bias, from the interaction 
between the particle and hole. We obtain the combined 
lb and eh FS correction from DFT calculations: 



SE. 



ih.lb 



(k) = 



; (fc)-A^ FT (fc), 



(4) 



where e g {k) = ej(k) — e^(fc) is the difference in band en- 



SEfj'(k) = AEjf A (k) - AEff~ LDA (k) 



(5) 



where the two excitation energies on the right are from 
LDA calculations paralleling the QMC, using the stan- 
dard and the FS functionals, respectively. The sum of 
Eqs. dU and J5]) gives the total FS correction SEij(k). 
The largest contribution is from <5£* h ' lb (/c), typically 
~ 0.10 eV at most k points in Si and diamond. Its largest 
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FIG. 2. (Color online) Many-body band structure of silicon 
from AFQMC (blue solid circles). GW £lj and LDA band 
structures are plotted by solid and dashed black lines, respec- 
tively. Available DMC results at high symmetry points, T, 
X, and L are indicated by green triangles; experimental 
values are given by red open circles. 



value in Si is 0.35 eV at the T point. In diamond, which 
has a large band gap, its largest value is 0.83 eV. The 2b 
correction is typically smaller; its largest value is 0.12 eV 
in Si and diamond, and 0.08 eV for the fundamental gap 
in ZnO. 

We then obtain a quasiparticle band structure e n (k) 
from a least-squares fit [lj| to the calculated many-body 
excitation energies AEij(k) = AE^ MC (k) + SE^k): 



^^(A^-(fc)-[ ej (fc)- ei (fc)]f 



(6) 



ergies from a standard DFT calculation, using a dense results 
fc-point grid, while AE® FT (k) is from self-consistent 
DFT calculations at k paralleling the QMC calculations 
in Eq. ([3]). We correct the 2b FS error, which along 
with the lb effect are substantially reduced here because 
A£<f c (fc) is an energy difference between two states 
within the same simulation cell, using an LDA functional 
specially parametrized lij 25| for FS calculations 



where i and j run over the occupied (v) and unoccupied 
(c) states, respectively. The highest occupied quasiparti- 
cle energy is set equal to the corresponding DFT eigen- 
value e m (k) = e m (k). 

The full many-body Si band structure, with FS correc- 
tion included, is shown in Fig. [2] and compared to LDA 
and GW, as well as to available DMC and experimental 

iEE|. 



Calculations for fee Si were done at the 
experimental lattice constant of 5.431 A. Both AFQMC 
and DFT calculations used a norm-conserving Kleinman- 
Bylander type separable nonlocal LDA pseudopotentials 
generated by the OPIUM code [27j . with a planewave 
cutoff E cut = 25 Ry. Using a |<0*) from LDA, AFQMC 
results correct the band gap problem and are in good 
agreement with experiment and with GW. The lowest 
band, which corresponds to the highest excitation ener- 
gies, tends to be ~ 1.5 eV too low. For an imaginary-time 
projection method, its quality can be expected to decline 
for higher excited states. Also the simple singlet \4>t) or 
the orthogonalization constraint may not be sufficient, as 
there are many states with similar energies. 

The many body band structure of diamond, with FS 
correction included, is given in Fig. EH tog ether with avail- 
able DMC and experimental results [16|, [21 1 . The lattice 
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FIG. 3. (Color online) Many-body band structure of diamond 
from AFQMC (blue solid circles). GW ^ and GGA band 
structures are plotted by solid and dashed black lines, respec- 
tively. Available DMC results at high symmetry points, L, 
X, and L [l6|], are indicated by green triangles; experimental 
values are given by red open circles. 



constant of 3.567 A was used. The calculations are sim- 
ilar to those in Si, except for a higher planewave cutoff 
^cut = 50 Ry and the use of GGA pseudopotential and 
trial wave function. We have verified that the calcula- 
tions are insensitive to the difference in at the level 
of LDA vs. GGA or a hybrid functional. The AFQMC 
results are again generally in very good agreement with 
GW and experiment. 

Accurately calculating the fundamental band gap of 
wurtzite structure ZnO is very challenging. DFT LDA 
and GGA underestimate the gap by almost 3eV. Even 
the generally more accurate GW method underestimates 
the gap by 0.8- 1.1 eV. There has been considerable dis- 
cussion about the importance of various convergence is- 
sues, choice of pseudopotentials, and additional approx- 
imations such as the plasmon-pole model and the inclu- 
sion of self-consistency in the GW calculation Re- 
sults are also sensitive to the choice of pseudopotential, 
since there is large spatial overlap among 3s, 3p and 3d 
wave functions of atomic Zn [9(. To properly treat this, 
our Zn GGA pseudopotential was constructed with a Ne- 
core, thereby fully correlating the semicore 3s, 3p along 
with the 3d states in the many-body calculation. Very 
conservative radial cutoffs of 0.83, 1.02, and 1.13 bohr 
were used for s, p, and d channels, respectively, result- 
ing in a large planewave cutoff energy of E cut = 180 Ry. 
The pseudopotential gives GGA optimized lattice param- 
eters of a = 3.279 and c = 5.284 A, and bulk modulus of 
128.7 GPa, all in good agreement with published GGA 
results [29j |. 

Our AFQMC calculations were done at the above GGA 
optimized geometry. While hybrid functionals seem to 
perform better in ZnO, we chose to use simple trial wave 
functions from GGA in our AFQMC calculations to avoid 
any parameter tuning. The singlet form |^) in Eq. @ 
was used for the excited state. The main calculations 



TABLE I. The calculated AFQMC band gap (eV) of wurtzite 
ZnO, compared to experiments [31-34]. Also shown are re- 
sults from GGA, hybrid DFT [lalirand GIL" [1, 11 [11 . 



GGA 


Hybrid 


GW 


AFQMC 


Expt. 


A gap 0.77 


3.32,2.90 


2.35,2.83,2.56 


3.26(16) 


3.30-3.57 



were done with a time-step of At = O.OlRy -1 . The Trot- 
ter error was then corrected by extrapolation from sepa- 
rate runs using a single-determinant \ip£). The calculated 
raw band gap was 2.54(14) eV. Since DFT severely under- 
estimates the band gap, the FS correction is less straight- 
forward in ZnO. At k = L, for example, the single /c-point 
self-consistent GGA calculation would yield a metallic 
ground state. We extrapolated SE^ ' (k) for k within 
0.1 of the L point, along two high symmetry lines. In do- 
ing so, we assume that the electron-hole effect does not 
introduce a discontinuity in the band dispersion, which 
is reasonable since it mainly relates to the simulation 



cell size. This yielded SE^ ch {0) = 0.58(08) eV. (Had we 

used the LDA, which has a smaller gap, <5_B I ? J b ' eh (0) would 
be about 0.41(04) eV.) Ths 2b correction is well-behaved: 
5Eff{0) = 0.08 eV. Adding the FS corrections yields our 
calculated band gap of 3.20(16) eV. 

The experimental equilibrium geometry is at a = 3.250 
and c = 5.207 A [3(| ■ To compare our band gap (for 
the GGA equilibrium geometry) to experiment, we apply 
a correction of +0.06 eV, which is the excitation energy 
difference given by GGA for the experimental and GGA- 
optimized geometries. (Hybrid B3LYP calculations give 
a correction of +0.10eV.) Table U compares our result 
with experimental values (3.30 (31 1 , 3.44 (32j, and 3.57 



eV 
381 



33, 34 



and those from recent calculations [8J, |3£ 
We note that, even with the small-core pseudopo- 
tential, there can be pseudoization and core-relaxation 
errors [391 ] . The precise effect cannot be determined with- 
out further many-body calculations. However, recent 
studies [35|, (3(| comparing all-electron and pseudopoten- 
tial GW calculations indicate that the effect is approxi- 
mately +0.27 eV for a pseudopotential of similar quality 
to the one we have adopted. This would increase the 
AFQMC result for the fundamental band gap in ZnO to 
~ 3.53(16) eV, within the range of experimental measure- 
ments [3lM34| listed in Table |U 

Various further improvements can be explored. We 
have used a planewave basis and norm-conserving pseu- 
dopotentials in these calculations. Other single-particle 
basis sets are straightforward to use, for example, with 
localized or natural orbitals. One could also work with a 
truncated set of orbitals from a lower-level of theory, or 
use a down-folded Hamiltonian directly to improve com- 
putational efficiency. The constraining virtual orbitals 
have been fixed in our calculations, but could potentially 
be allowed to dynamically evolve in some way in the 
propagation. It is reassuring that DFT trial wave func- 
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tions have worked well. Clearly more elaborate multi- 
determinant wave functions can be used, for example, 
from diagonalization in a subspace formed by single and 
double excitations to conduction orbitals. 

In summary, we have presented an AFQMC approach 
for the calculation of electronic excitations in solids, 
introducing an orbitally-based orthogonalization con- 
straint in the phaseless AFQMC framework to stabi- 
lize the projection of excited states. Simple trial wave 
functions directly from DFT calculations were used for 
the constraint. Detailed many-body quasiparticle band 
structures can be calculated. In prototypical semiconduc- 
tors (Si and diamond) , the calculated band structures are 
in good agreement with those from GW calculations and 
with experiment. In the more strongly correlated and 
challenging wurtzite ZnO crystal, the calculated funda- 
mental gap is in excellent agreement with the latest ex- 
perimental data. The method is non-perturbative and 
free of empirical parameters, offering a possible path for 
general computations in correlated materials. 
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